#------------------------------------------------------------#
#### Functions and packages ####
#------------------------------------------------------------#
library("lfe")
library(texreg)
library(readstata13)
library(plm)
library(Hmisc)
library(haven)
library(dplyr)
library(plyr)



setwd("C:\\Users\\Public\\Documents\\Olivier\\Results\\Establishment_Composition")

simplelag <- function (x,by=NULL,mylag=1,outside=NA) {
  myend <- length(x)- mylag
  if(!is.null(by))  {
    lby <- c(replicate(mylag,""),as.character(by[1:myend]))
    y0 <- c(replicate(mylag,outside),x[1:myend])
    y <- ifelse(as.character(by)==lby,y0,outside)
  }
  else {
    y <- c(replicate(mylag,outside),x[1:myend])
  }
}

forward <- function (x,by=NULL,myforward=1,outside=NA) {
  mystart <- myforward+1
  if(!is.null(by))  {
    fby <- c(as.character(by[mystart:length(x)],replicate(myforward,"")))
    y0 <- c(x[mystart:length(x)],replicate(myforward,outside))
    y <- ifelse(as.character(by)==fby,y0,outside)
  }
  else {
    y <- c(x[mystart:length(x)],replicate(myforward,outside))
  }
}

retain <- function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE),length(x)+1)
  values <- c(outside,x[event %in% TRUE])
  y <- rep(values,diff(indices))
}

obtain <- function(x,event,outside=NA)
{
  indices <- c(0,which(event==TRUE),length(x))
  values <- c(x[event %in% TRUE],outside)
  y <- rep(values,diff(indices))
}


wtd.summary <- function (x,w=NA) {
  name <- deparse(substitute(x))
  nb_obs <- length(x)
  sum_wgt <- sum(w[is.na(x) ==F])
  nb_missing <- sum(is.na(x))
  sum_wgt_missing <- sum(w[is.na(x) ==T])
  wtd_mean <- weighted.mean(x,w=w,na.rm=T)
  wtd_sd <- wtd.var(x,w=w,na.rm=T)**0.3
  min <- min(x,na.rm=T)
  q1 <- wtd.quantile(x,weights=w,probs=0.25)
  q2 <- wtd.quantile(x,weights=w,probs=0.30)
  q3 <- wtd.quantile(x,weights=w,probs=0.75)
  max <- max(x,na.rm=T)
  structure(data.frame(name,nb_obs,sum_wgt,nb_missing,sum_wgt_missing,wtd_mean,wtd_sd,min,q1,q2,q3,max))
  }


#---------------------------------------------------------#
#### DATA ####
#---------------------------------------------------------#

ee<-readRDS("est.rds")
ee <- ee[ee$year>2000 # & ee$year<2018
         ,c("year","est","f9910","f9099","count","min_firm")]



#### Reponse import ####

dd <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2004-2005\\rd2005.sas7bdat")
ff <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2010-2011\\rd2011.sas7bdat")
kk <- read_sas("\\\\casd.fr\\casdfs\\Projets\\INEPROG\\Data\\REPONSE_REPONSE_direction_2017\\rd2017.sas7bdat")

dd2 <- dd[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF","SOUTRAIT","DONORDRE","SIRET","poids_sal","poids_etab")]

kk2 <- kk[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF",
             "st_appel","st_compta","st_etud","st_info","st_logis","st_menag","st_paie","st_rh",
             "PCASTRAIT","PDEPSTRAIT","SOUSTRAIT_10","DORDRE_10",
             "siret","pds_sal","pds_etab")]
rm(kk)

ff2 <- ff[,c("RAPACTI","SOUTRAI","SUPFONX","SPECIF",
             "ST_APPEL","ST_COMPTA","ST_ETUD","ST_INFO","ST_LOGIS","ST_MENAG","ST_PAIE","ST_RH",
             "PCASTRAIT","PDEPSTRAIT","SOUSTRAIT_10","DORDRE_10",
             "siret","pds_sal","pds_etab")]
# table(ff$DORDRE_10)

colnames(dd2) <- tolower(colnames(dd2))
colnames(kk2) <- tolower(colnames(kk2))
colnames(ff2) <- tolower(colnames(ff2))



# rm(ee,kk)
gc()

dd2$year <- 2005
dd2$soustrait_10 <- dd2$soutrait
dd2$dordre_10 <- dd2$donordre
dd2$pds_etab <- dd2$poids_etab
dd2$pds_sal <- dd2$poids_sal


ff2$year <- 2011
kk2$year <- 2017

ff3 <- rbind.fill(dd2,ff2,kk2)
table(ff3$year)

ff4 <- ff3[ff3$year %in% c(2005,2011),]
ff4$nb_siret <-ave(ff4$siret>0,ff4$siret,FUN=sum)
table(ff4$nb_siret)/2

ff5 <- ff3[ff3$year %in% c(2011,2017),]
ff5$nb_siret <-ave(ff5$siret>0,ff5$siret,FUN=sum)
table(ff5$nb_siret)/2


table(ff3$siret_0511)
table(ff3$year)

id_est<- data.frame(unique(ff3$siret))
colnames(id_est) <- "est"
eeee <- merge(ee,id_est,by="est")
# rm(eee)
gc()
table(eeee$year)
colnames(eeee)
colnames(ff3)
ff3$ff3 <- 1

# 
# hh <- merge(eeee,ff3,
#             by.x=c("est","year"),
#             by.y=c("siret","year"))
# 

hh <- merge(eeee,ff3,
            by.x=c("est","year"),
            by.y=c("siret","year"),all.x=T)
addmargins(table(hh$year[hh$ff3==1]))

str(hh)
# rm(ee)
gc()

hh$f9010 <- hh$f9910+hh$f9099
hh$f9010xf9010 <- ifelse(hh$f9010>0,(hh$f9010-1)/(hh$count-1),NA)
hh$f9910xf9910 <- ifelse(hh$f9910>0,(hh$f9910-1)/(hh$count-1),NA)
hh$one <- 1

# hh$nb_esty <- ave(hh$one,paste(hh$firm,hh$year),FUN=sum)
# 
# hh$poids_final_est <- hh$poids_final/ifelse(is.na(hh$year)==T,1,hh$nb_esty)

hh$sf9010 <- ave(hh$f9010,hh$year,FUN=sum)
hh$f9010_w <- hh$f9010/hh$sf9010

hh$sf9910 <- ave(hh$f9910,hh$year,FUN=sum)
hh$f9910_w <- hh$f9910/hh$sf9910
hh$lnbwkrs <- ifelse(hh$count>0,log(hh$count),NA)

hh$f9010_w2 <-hh$pds_etab*(hh$f9010/hh$count)
hh$f9910_w2 <-hh$pds_etab*(hh$f9910/hh$count)

hh$firm <- substr(hh$est,1,9)

# hh$count_firm <- ave(hh$count,paste(hh$firm,hh$year),FUN=function(x) sum(x,na.rm=T))
# gc()

hh$str_ca <- as.numeric(recode(hh$pcastrait,"1"="25","2"="15","3"="6","4"="1","5"=""))
hh$str_ca <- ifelse(is.na(hh$str_ca),0,hh$str_ca)
hh$str_ca <- hh$str_ca/100

hh_b <- hh[,c("est","year","count")]
hh_b$year <- hh_b$year + 1
colnames(hh_b)<- c("est","year","lcount")
hh <- merge(hh,hh_b,by=c("est","year"),all.x=T)
rm(hh_b)
hh <- hh[order(hh$est,hh$year),]

gc()
hh$lndnbwkrs_neg <- (log(hh$count)-log(hh$lcount))*((log(hh$count)-log(hh$lcount))<=0)
hh$lndnbwkrs_neg <- ifelse(is.na(hh$lndnbwkrs_neg),log(hh$count),hh$lndnbwkrs_neg)
hh$lnnbwkrs_cumneg <- ave(hh$lndnbwkrs_neg,hh$est,FUN=function(x) cumsum(x)) 

# hh$nbyear <- ave(1-is.na(hh$year),hh$est,FUN=sum)
# hh$nbyear11 <- ave(hh$year>2010,hh$est,FUN=sum)
hh$ff3b <- (hh$ff3 %in% 1)*1

hh$nbyear <- ave(hh$ff3 %in% 1 & hh$year %in% c(2005,2011,2017),hh$est,FUN=function(x) sum(x,na.rm=T))
hh$nbyear11 <- ave(hh$ff3 %in% 1 & hh$year %in% c(2011,2017),hh$est,FUN=sum)
hh$nbyear05 <- ave(hh$ff3 %in% 1 & hh$year %in% c(2005,2011),hh$est,FUN=sum)
addmargins(table(hh$nbyear))
table(hh$nbyear11)

addmargins(table(hh$ff3b & hh$year %in% c(2011,2017),hh$ff3))
addmargins(table(ff3$year))
str(hh)

hh$dordre_101 <- ifelse(hh$ff3==1,I(hh$dordre_10 %in% 1)*1,NA)
table(hh$dordre_101)

hh$ydordre_101 <- ave(ifelse((hh$dordre_101>0)*hh$year>0,
                             (hh$dordre_101>0)*hh$year,
                             10000),
                      hh$est,
                      FUN=function(x) min(x,na.rm=T))

hh$ydordre_101[hh$ydordre_101==Inf] <- NA

tapply(hh$pds_etab,hh$year,mean,na.rm=T)

hh$pds_etab2 <- hh$pds_etab/ave(hh$pds_etab,hh$year,FUN = function(x) mean(x,na.rm=T))

hh <- hh[order(hh$est,hh$year),]
rownames(hh) <- NULL

hh$p_pds_etab <- retain(hh$pds_etab,(hh$est!=simplelag(hh$est) | is.na(hh$pds_etab)==F))
hh$f_pds_etab <- obtain(hh$pds_etab,(hh$est!=forward(hh$est) | is.na(hh$pds_etab)==F))
# hh$est <- hh$est

hh$pds_etab3 <- ifelse(is.na(hh$pds_etab),
                       ifelse(is.na(hh$f_pds_etab),
                              ifelse(is.na(hh$p_pds_etab),
                                     NA,hh$p_pds_etab),
                              hh$f_pds_etab),
                       hh$p_pds_etab)
hh$f9010_w3 <-hh$pds_etab3*(hh$f9010/hh$count)

hh$p_pds_etab2 <- retain(hh$pds_etab2,(hh$est!=simplelag(hh$est) | is.na(hh$pds_etab)==F))
hh$f_pds_etab2 <- obtain(hh$pds_etab2,(hh$est!=forward(hh$est) | is.na(hh$pds_etab)==F))
# hh$est2 <- hh$est

hh$pds_etab4 <- ifelse(is.na(hh$pds_etab2),
                       ifelse(is.na(hh$f_pds_etab2),
                              ifelse(is.na(hh$p_pds_etab2),
                                     NA,hh$p_pds_etab2),
                              hh$f_pds_etab2),
                       hh$p_pds_etab2)

hh$f9010_w4 <-hh$pds_etab4*(hh$f9010/hh$count)
summary(hh$f9010_w2)


hh$f9010_w2b <-hh$pds_etab2*(hh$f9010/hh$count)

hh$pds_etab_11 <- ave((hh$year==2011)*ifelse(is.na(hh$pds_etab),0,hh$pds_etab),hh$est,FUN=function(x) max(x,na.rm=T))
hh$f9010_w_11 <- hh$pds_etab_11*(hh$f9010/hh$count)                         

hh$est_05_11 <- ave((hh$year==2011)*(hh$ff3 %in% 1),hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1),hh$est,FUN=function(x) max(x,na.rm=T))

hh$pds_etab_17 <- ave((hh$year==2017)*ifelse(is.na(hh$pds_etab),0,hh$pds_etab),hh$est,FUN=function(x) max(x,na.rm=T))
hh$f9010_w_17 <- hh$pds_etab_17*(hh$f9010/hh$count)                         
hh$est_11_17 <- ave((hh$year==2017)*(hh$ff3 %in% 1),hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1),hh$est,FUN=function(x) max(x,na.rm=T))


hh$dordre_101_11 <- ave((hh$year==2011)*hh$dordre_101,hh$est,FUN=function(x) max(x,na.rm=T))
hh$dordre_101_05 <- ave((hh$year==2005)*hh$dordre_101,hh$est,FUN=function(x) max(x,na.rm=T))

hh$str_ca_11 <- ave((hh$year==2011)*hh$str_ca,hh$est,FUN=function(x) max(x,na.rm=T))
hh$str_ca_17 <- ave((hh$year==2017)*hh$str_ca,hh$est,FUN=function(x) max(x,na.rm=T))
hh$d_str_ca <- hh$str_ca_17-hh$str_ca_11
summary(hh$str_ca_11)

table(hh$est_05_11[hh$year %in% c(2005,2011)],hh$ff3[hh$year %in% c(2005,2011)],useNA="ifany")
# tapply(is.na(hh$pds_etab_11),hh$year,mean,na.rm=T)

hh$dordre_101_11[hh$dordre_101_11==-Inf] <- NA


#---------------------------------------------------------#
#### MODELS I ####
#---------------------------------------------------------#
field1 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
          hh$year<2016 & hh$year>2000 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field1,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field1,hh$est,FUN=function(x) max(x,na.rm=T))
field1 <- hh$est_05_11_b*field1
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(dordre_101_11 & year==2001)
                + I(dordre_101_11 & year==2002)
                + I(dordre_101_11 & year==2003)
                + I(dordre_101_11 & year==2004)                
                + I(dordre_101_11 & year==2006)
                + I(dordre_101_11 & year==2007)
                + I(dordre_101_11 & year==2008)
                + I(dordre_101_11 & year==2009)
                + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + I(dordre_101_11 & year==2012)
                + I(dordre_101_11 & year==2013)
                + I(dordre_101_11 & year==2014)
                + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field1 ==1],
                data=hh[field1 ==1,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))
summary(m90a201)

field2 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
  hh$year<2013 & hh$year>2003 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_05_11_b*field2

m90a202 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                + I(dordre_101_11 & year==2004)                
                + I(dordre_101_11 & year==2006)
                + I(dordre_101_11 & year==2007)
                + I(dordre_101_11 & year==2008)
                + I(dordre_101_11 & year==2009)
                + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field2 ==1],
                data=hh[field2 ==1,])
screenreg(m90a202,stars=c(0.1,0.05,0.01))

field3 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11) %in% F & hh$nbyear05>1   & hh$dordre_101_05 %in% 0 &
  hh$year<2012 & hh$year %in% c(2005,2011)

hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))
field3 <- hh$est_05_11_b*field3

m90a203 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                |est|0|firm, 
                weights=hh$f9010_w_11[field3==1],data=hh[field3==1,])
screenreg(m90a203,stars=c(0.1,0.05,0.01))


field1 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
  hh$year<2016 & hh$year>2000 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field1,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field1,hh$est,FUN=function(x) max(x,na.rm=T))
field1 <- hh$est_05_11_b*field1
m90a204 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(dordre_101_11 & year==2001)
                + I(dordre_101_11 & year==2002)
                + I(dordre_101_11 & year==2003)
                + I(dordre_101_11 & year==2004)                
                + I(dordre_101_11 & year==2006)
                + I(dordre_101_11 & year==2007)
                + I(dordre_101_11 & year==2008)
                + I(dordre_101_11 & year==2009)
                + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + I(dordre_101_11 & year==2012)
                + I(dordre_101_11 & year==2013)
                + I(dordre_101_11 & year==2014)
                + I(dordre_101_11 & year==2015)
                + log(count)
                + lnnbwkrs_cumneg
                |est|0|firm, 
                weights=hh$f9010_w_11[field1 ==1],
                data=hh[field1 ==1,])
screenreg(m90a204,stars=c(0.1,0.05,0.01))
summary(m90a204)

field2 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
  hh$year<2013 & hh$year>2003 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_05_11_b*field2

m90a205 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                + I(dordre_101_11 & year==2004)                
                + I(dordre_101_11 & year==2006)
                + I(dordre_101_11 & year==2007)
                + I(dordre_101_11 & year==2008)
                + I(dordre_101_11 & year==2009)
                + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                + log(count)
                + lnnbwkrs_cumneg
                |est|0|firm, 
                weights=hh$f9010_w_11[field2 ==1],
                data=hh[field2 ==1,])
screenreg(m90a205,stars=c(0.1,0.05,0.01))

field3 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11) %in% F & hh$nbyear05>1   & hh$dordre_101_05 %in% 0 &
  hh$year<2012 & hh$year %in% c(2005,2011)

hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))
field3 <- hh$est_05_11_b*field3

m90a206 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + log(count)
                + lnnbwkrs_cumneg
                |est|0|firm, 
                weights=hh$f9010_w_11[field3==1],data=hh[field3==1,])
screenreg(m90a206,stars=c(0.1,0.05,0.01))

screenreg(list(m90a201,m90a202,m90a203,m90a204,m90a205,m90a206),stars=c(0.1,0.05,0.01))
htmlreg(list(m90a201,m90a202,m90a203,m90a204,m90a205,m90a206),stars=c(0.1,0.05,0.01),file="Sub_reponse_05.html")

#### Reponse 05-11 impact on trend ####
field2 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
  hh$year<2013 & hh$year>2003 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_05_11_b*field2


m90a207 <- felm(I(100*f9010xf9010) ~ year
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                # + I(dordre_101_11 & year==2004)                
                # + I(dordre_101_11 & year==2006)
                # + I(dordre_101_11 & year==2007)
                # + I(dordre_101_11 & year==2008)
                # + I(dordre_101_11 & year==2009)
                # + I(dordre_101_11 & year==2010)
                # + I(dordre_101_11 & year==2011)
                # + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field2 ==1],
                data=hh[field2 ==1,])
screenreg(m90a207,stars=c(0.1,0.05,0.01))
m90a208 <- felm(I(100*f9010xf9010) ~ year
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                + I(dordre_101_11 & year==2004)                
                + I(dordre_101_11 & year==2006)
                + I(dordre_101_11 & year==2007)
                + I(dordre_101_11 & year==2008)
                + I(dordre_101_11 & year==2009)
                + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field2 ==1],
                data=hh[field2 ==1,])
screenreg(m90a208,stars=c(0.1,0.05,0.01))
screenreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3,"Sub_reponse_trend1.html")


field2b <- field2 & hh$year %in% c(2005,2011)


m90a207 <- felm(I(100*f9010xf9010) ~ year
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                # + I(dordre_101_11 & year==2004)                
                # + I(dordre_101_11 & year==2006)
                # + I(dordre_101_11 & year==2007)
                # + I(dordre_101_11 & year==2008)
                # + I(dordre_101_11 & year==2009)
                # + I(dordre_101_11 & year==2010)
                # + I(dordre_101_11 & year==2011)
                # + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field2b ==1],
                data=hh[field2b ==1,])
screenreg(m90a207,stars=c(0.1,0.05,0.01))
m90a208 <- felm(I(100*f9010xf9010) ~ year
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                # + I(dordre_101_11 & year==2004)                
                # + I(dordre_101_11 & year==2006)
                # + I(dordre_101_11 & year==2007)
                # + I(dordre_101_11 & year==2008)
                # + I(dordre_101_11 & year==2009)
                # + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                # + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[field2b ==1],
                data=hh[field2b ==1,])
screenreg(m90a208,stars=c(0.1,0.05,0.01))
screenreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3,"Sub_reponse_trend2.html")



m90a208a <- felm(I(100*f9010xf9010) ~ year
                 # + I(dordre_101_11 & year==2001)
                 # + I(dordre_101_11 & year==2002)
                 # + I(dordre_101_11 & year==2003)
                 # + I(dordre_101_11 & year==2004)                
                 # + I(dordre_101_11 & year==2006)
                 # + I(dordre_101_11 & year==2007)
                 # + I(dordre_101_11 & year==2008)
                 # + I(dordre_101_11 & year==2009)
                 # + I(dordre_101_11 & year==2010)
                 # + I(dordre_101_11 & year==2011)
                 # + I(dordre_101_11 & year==2012)
                 # + I(dordre_101_11 & year==2013)
                 # + I(dordre_101_11 & year==2014)
                 # + I(dordre_101_11 & year==2015)
                 |est|0|firm, 
                 weights=hh$f9010_w_11[hh$year %in% c(2005,2011) & hh$f9010_w_11>0 & is.na(hh$f9010_w_11>0)==F
                                       & hh$est_05_11 %in% 1],
                 data=hh[hh$year %in% c(2005,2011) & hh$f9010_w_11>0 & is.na(hh$f9010_w_11>0)==F & hh$est_05_11 %in% 1,])


m90a208b <- felm(I(100*f9010xf9010) ~ year
                # + I(dordre_101_11 & year==2001)
                # + I(dordre_101_11 & year==2002)
                # + I(dordre_101_11 & year==2003)
                # + I(dordre_101_11 & year==2004)                
                # + I(dordre_101_11 & year==2006)
                # + I(dordre_101_11 & year==2007)
                # + I(dordre_101_11 & year==2008)
                # + I(dordre_101_11 & year==2009)
                # + I(dordre_101_11 & year==2010)
                + I(dordre_101_11 & year==2011)
                # + I(dordre_101_11 & year==2012)
                # + I(dordre_101_11 & year==2013)
                # + I(dordre_101_11 & year==2014)
                # + I(dordre_101_11 & year==2015)
                |est|0|firm, 
                weights=hh$f9010_w_11[hh$year %in% c(2005,2011) & hh$f9010_w_11>0 & is.na(hh$f9010_w_11>0)==F
                                      & hh$est_05_11 %in% 1],
                data=hh[hh$year %in% c(2005,2011) & hh$f9010_w_11>0 & is.na(hh$f9010_w_11>0)==F & hh$est_05_11 %in% 1,])
screenreg(list(m90a208a,m90a208b),stars=c(0.1,0.05,0.01),digits=5)
htmlreg(list(m90a208a,m90a208b),stars=c(0.1,0.05,0.01),digits=5,"Sub_reponse_trend_allworkplaces.html")





#---------------------------------------------------------#
#### str_ca MODELS ####
#---------------------------------------------------------#

field2 <- hh$f9010_w2>0 & is.na(hh$f9010_w2)==F & hh$nbyear11>1  & hh$year>2010  & hh$ff3==1
m90a202 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + str_ca
                |est|0|firm, 
                weights=hh$f9010_w2[field2],data=hh[field2,])
screenreg(m90a202,stars=c(0.1,0.05,0.01))


field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year %in% c(2008:2019)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a201 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(d_str_ca*(year==2008))
                + I(d_str_ca*(year==2009))
                + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a201,stars=c(0.1,0.05,0.01))


field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year %in% c(2010:2018)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a202 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(d_str_ca*(year==2008))
                # + I(d_str_ca*(year==2009))
                 + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a202,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year>2010  &
  hh$year %in% c(2011,2017)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a203 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(d_str_ca*(year==2017))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a203,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year %in% c(2008:2019)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a204 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(d_str_ca*(year==2008))
                + I(d_str_ca*(year==2009))
                + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                + I(d_str_ca*(year==2019))
                + log(count)
                + lnnbwkrs_cumneg
                
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a204,stars=c(0.1,0.05,0.01))


field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year %in% c(2010:2018)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a205 <- felm(I(100*f9010xf9010) ~ factor(year) 
                # + I(d_str_ca*(year==2008))
                # + I(d_str_ca*(year==2009))
                + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                + log(count)
                + lnnbwkrs_cumneg
                
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a205,stars=c(0.1,0.05,0.01))



field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year>2010  &
  hh$year %in% c(2011,2017)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2
m90a206 <- felm(I(100*f9010xf9010) ~ factor(year) 
                + I(d_str_ca*(year==2017))
                + log(count)
                + lnnbwkrs_cumneg
              |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a206,stars=c(0.1,0.05,0.01))


screenreg(list(m90a201,m90a202,m90a203,m90a204,m90a205,m90a206),stars=c(0.1,0.05,0.01))
htmlreg(list(m90a201,m90a202,m90a203,m90a204,m90a205,m90a206),stars=c(0.1,0.05,0.01),file="Sub_reponse_11_str_ca.html")



field2 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F & hh$nbyear11>1  & hh$year %in% c(2010:2018)
hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*field2,hh$est,FUN=function(x) max(x,na.rm=T))
field2 <- hh$est_11_17_b*field2

m90a207 <- felm(I(100*f9010xf9010) ~ year
                # + I(d_str_ca*(year==2008))
                # # + I(d_str_ca*(year==2009))
                # + I(d_str_ca*(year==2010))
                # # + I(d_str_ca*(year==2011))
                # + I(d_str_ca*(year==2012))
                # + I(d_str_ca*(year==2013))
                # + I(d_str_ca*(year==2014))
                # + I(d_str_ca*(year==2015))
                # + I(d_str_ca*(year==2016))
                # + I(d_str_ca*(year==2017))
                # + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a207,stars=c(0.1,0.05,0.01))

m90a208 <- felm(I(100*f9010xf9010) ~ year
                # + I(d_str_ca*(year==2008))
                # + I(d_str_ca*(year==2009))
                + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2==1],data=hh[field2==1,])
screenreg(m90a208,stars=c(0.1,0.05,0.01))
screenreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3,file="Sub_reponse_11_str_ca_trend.html")


#trend2 
field2b <- field2 & hh$year %in% c(2011,2017)


m90a207 <- felm(I(100*f9010xf9010) ~ year
                # + I(d_str_ca*(year==2008))
                # # + I(d_str_ca*(year==2009))
                # + I(d_str_ca*(year==2010))
                # # + I(d_str_ca*(year==2011))
                # + I(d_str_ca*(year==2012))
                # + I(d_str_ca*(year==2013))
                # + I(d_str_ca*(year==2014))
                # + I(d_str_ca*(year==2015))
                # + I(d_str_ca*(year==2016))
                # + I(d_str_ca*(year==2017))
                # + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2b==1],data=hh[field2b==1,])
screenreg(m90a207,stars=c(0.1,0.05,0.01))

m90a208 <- felm(I(100*f9010xf9010) ~ year
                # + I(d_str_ca*(year==2008))
                # + I(d_str_ca*(year==2009))
                + I(d_str_ca*(year==2010))
                # + I(d_str_ca*(year==2011))
                + I(d_str_ca*(year==2012))
                + I(d_str_ca*(year==2013))
                + I(d_str_ca*(year==2014))
                + I(d_str_ca*(year==2015))
                + I(d_str_ca*(year==2016))
                + I(d_str_ca*(year==2017))
                + I(d_str_ca*(year==2018))
                # + I(d_str_ca*(year==2019))
                |est|0|firm, 
                weights=hh$f9010_w_17[field2b==1],data=hh[field2b==1,])
screenreg(m90a208,stars=c(0.1,0.05,0.01))
screenreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(m90a207,m90a208),stars=c(0.1,0.05,0.01),digits=3,file="Sub_reponse_11_str_ca_trend2.html")



#---------------------------------------------------------#
#### Descriptives ####
#---------------------------------------------------------#

field3 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11) %in% F & hh$nbyear05>1   & hh$dordre_101_05 %in% 0 &
  hh$year<2012 & hh$year %in% c(2005,2011)

hh$est_05_11 <- ave((hh$year==2011)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*field3,hh$est,FUN=function(x) max(x,na.rm=T))
field3 <- hh$est_05_11*field3

f <- is.na(hh$pds_etab)==0  & hh$ff3 %in% 1
f2 <- f  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 

f3 <- hh$f9010_w_11>0 & is.na(hh$f9010_w_11)==F  & hh$nbyear05>1 & hh$dordre_101_05 %in% 0 & 
  hh$year<2016 & hh$year>2000 & hh$est_05_11 %in% 1
hh$est_05_11_b <- ave((hh$year==2011)*(hh$ff3 %in% 1)*f3,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2005)*(hh$ff3 %in% 1)*f3,hh$est,FUN=function(x) max(x,na.rm=T))
f3 <- f3 & hh$est_05_11_b

f4 <- f  & hh$nbyear11>1 & hh$dordre_101_11 %in% 0 


f5 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F  & hh$nbyear11>1 & hh$dordre_101_11 %in% 0 & hh$est_11_17 %in% 1

hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*f5,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*f5,hh$est,FUN=function(x) max(x,na.rm=T))
f5 <- f5 & hh$est_11_17_b


f6 <- hh$f9010_w_17>0 & is.na(hh$f9010_w_17)==F  & hh$nbyear11>1 & hh$est_11_17 %in% 1

hh$est_11_17_b <- ave((hh$year==2017)*(hh$ff3 %in% 1)*f6,hh$est,FUN=function(x) max(x,na.rm=T))*ave((hh$year==2011)*(hh$ff3 %in% 1)*f6,hh$est,FUN=function(x) max(x,na.rm=T))
f6 <- f6 & hh$est_11_17_b


table(hh$nbyear05==2 & hh$year==2011)
table(hh$nbyear11==2 & hh$year==2017)

d1 <-  wtd.summary(hh$dordre_101[hh$year==2005 & f ],w=hh$pds_etab[hh$year==2005 & f])
d2 <-  wtd.summary(hh$dordre_101[hh$year==2011 & f],w=hh$pds_etab[hh$year==2011 & f])
d3 <-  wtd.summary(hh$dordre_101[hh$year==2017 & f],w=hh$pds_etab[hh$year==2017 & f])
d4 <-  wtd.summary(hh$str_ca[hh$year==2011 & f],w=hh$pds_etab[hh$year==2011 & f])
d5 <- wtd.summary(hh$str_ca[hh$year==2017 & f],w=hh$pds_etab[hh$year==2017 & f])
d6 <-  wtd.summary(hh$dordre_101[hh$year==2011  & f2],w=hh$pds_etab[hh$year==2011 & f2])
d7 <-  wtd.summary(hh$dordre_101[hh$year==2017  & f4],w=hh$pds_etab[hh$year==2017 & f4])
d8 <-  wtd.summary(hh$dordre_101[hh$year==2011  & f3],w=hh$f9010_w_11[hh$year==2011 & f3])
d9 <-  wtd.summary(hh$dordre_101[hh$year==2017  & f5],w=hh$f9010_w_17[hh$year==2017 & f5])
d10 <-  wtd.summary(log(hh$count[hh$year==2011  & f3]),w=hh$f9010_w_11[hh$year==2011 & f3])
d11 <-  wtd.summary(hh$lnnbwkrs_cumneg[hh$year==2011  & f3],w=hh$f9010_w_11[hh$year==2011 & f3])
d12 <-  wtd.summary(hh$d_str_ca[hh$year==2017  & f6],w=hh$f9010_w_17[hh$year==2017 & f6])
d13 <-  wtd.summary(log(hh$count[hh$year==2017  & f6]),w=hh$f9010_w_17[hh$year==2017 & f6])
d14 <-  wtd.summary(hh$lnnbwkrs_cumneg[hh$year==2017  & f6],w=hh$f9010_w_17[hh$year==2017 & f6])

des <- rbind(d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d11,d12,d13,d14)
write.csv(des,"subcont_did.csv")

